---
title: Appendix - PureSpectrum Survey
author: Weifang Xu, Taylor Chewning, and Qing Wang
date: October 8, 2024
output: pdf_document
fontsize: 11 pt
header-includes:
  \usepackage[T1]{fontenc}
  \usepackage[utf8]{inputenc}
  \usepackage{newpxtext,newpxmath}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
## Clean the working environment and set up the working directory 
rm(list = ls())
setwd("/Users/qingwang/Downloads/Data Replication")

# load the libraries
library(readr)
library(tidyverse)
library(kableExtra)
library(haven)
library(ggthemes)
library(mediation)
library(flexplot)
library(sandwich)
library(lmtest)
library(texreg)
library(boot)
library(xtable)
library(modelsummary)
library(marginaleffects)
library(ggpubr)
library(pBrackets)
library(lemon)
library(arm)
library(gplots)
library(extrafont)
library(broom)
library(multcomp)
library(ggeffects)
# install the compactr package
# url <- "https://cran.r-project.org/src/contrib/Archive/compactr/compactr_0.1.tar.gz"
# install.packages(url, repos = NULL, type = "source")
library(compactr)

# load the cleaned dataset
df_clean <- read_rds("PureSpectrum/clean_data_PureSpectrum.rds")
```


```{r}
#### Table S1: Sample Demographics in Comparison with Census Benchmarks (PureSpectrum Survey) ####

# benchmark demographic data is obtained from Table S0101 of the 2021 American Community Survey
# link to the survey: https://data.census.gov/table/ACSST1Y2021.S0101?q=S0101

# calculate the demographic of the PureSpectrum sample
# sex
male_percent <- df_clean %>%
  group_by(male) %>%
  summarise(percentage = round(n() / nrow(df_clean) * 100, 1))
male_percent

# age
# age ==1, 18–29
# age ==2, 30–39
# age ==3, 40–49
# age ==4, 50–59
# age ==5, 60–69
# age ==6, 70 or above

age_percent <- df_clean %>%
  group_by(age) %>%
  summarise(percentage = round(n() / nrow(df_clean) * 100, 1))
age_percent

# race
# race ==1, White
# race ==2, Black or African American
# race ==3, Asian American
# race ==4, Hispanic
# race ==5, Native American
# race ==6, Other

race_percent <- df_clean %>%
  group_by(race) %>%
  summarise(percentage = round(n() / nrow(df_clean) * 100, 1))
race_percent

# calculate the "Other" race category
race_other = race_percent[3,2] + race_percent[5,2] + race_percent[6,2]
race_other

```


```{r}
#### Table S2: Mean of Dependent Variable by Treatment Condition (PureSpectrum Survey) ####

# excluding the obs with NA in attack
df_clean_complete <- df_clean %>% filter(complete.cases(attack))

mean_dv <- df_clean_complete %>% 
  group_by(exp_4) %>%
  summarise(mean = mean(attack, na.rm = TRUE),
            sd = sd(attack, na.rm = TRUE),
            n = n()) %>%
  mutate_if(is.numeric, round, digits=2)

mean_dv %>%
  kbl(caption = "Mean of Dependent Variable by Treatment Condition", # Adding caption  
      format = "latex") # Output format = latex

```


```{r}
#### Figure S1: Conditional Average Treatment Effects (CATE) by Party ####

# Create 4-category treatment variable  
df_heter <- df_clean %>% 
  mutate(treatment_cat = ifelse(hmrts == 1 & alliance == 0, 1,
                              ifelse(hmrts == 1 & alliance == 1, 2,
                                     ifelse(hmrts == 0 & alliance == 1, 3,
                                            ifelse(hmrts == 0 & alliance == 0, 4, NA)))),
         treatment_cat = factor(treatment_cat, levels = c("1", "2", "3", "4"),
                              labels = c("Non-Ally, Violates HR", "Ally, Violates HR", "Ally, Protects HR", "Non-Ally, Protects HR")))%>%       
  glimpse()

#Run Interaction Model; Filter out Independents
mod <- lm(attack ~ treatment_cat*party + male + age_cat + edu4 + inc + white, data = subset(df_heter, party != "Independent"))
summary(mod)

#Store regression estimates in df
CATE_summary_df <- broom::tidy(mod) %>%
  #Filter out irrelevant terms  
  filter(grepl('treatment_cat', term))%>%
  #Rename terms  
  mutate(term = factor(term,
                       levels = c("treatment_catAlly, Violates HR", 
                                  "treatment_catAlly, Protects HR", 
                                  "treatment_catNon-Ally, Protects HR", 
                                  "treatment_catAlly, Violates HR:partyRepublican",
                                  "treatment_catAlly, Protects HR:partyRepublican",
                                  "treatment_catNon-Ally, Protects HR:partyRepublican"),
                       labels = c("CATE on Democrats", #Effect of Treatment on DV when Republican == 0
                                  "CATE on Democrats", 
                                  "CATE on Democrats",
                                  "Difference in CATEs", 
                                  "Difference in CATEs", 
                                  "Difference in CATEs")))%>%
  mutate(treatment_cat = c(1,2,3,1,2,3),
         treatment_cat = factor(treatment_cat,
                                levels = 1:3,
                                labels = c("Ally,\nViolates HR",
                                           "Ally,\nProtects HR",
                                           "Non-Ally,\nProtects HR")))%>%
  glimpse()



#### Calculate CATE for Republicans 

temp_estimate_1 <- glht(mod, linfct = c("`treatment_catAlly, Violates HR` + `treatment_catAlly, Violates HR:partyRepublican` = 0"))
temp_df_1 <- 
  as.data.frame(matrix(c("CATE on Republicans", 
                         confint(temp_estimate_1)$confint[ , c("Estimate", "lwr", "upr")], 
                         "attack", "Ally,\nViolates HR"), 
                       nrow = 1))
temp_df_1 <- temp_df_1 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)


temp_estimate_2 <- glht(mod, linfct = c("`treatment_catAlly, Protects HR` + `treatment_catAlly, Protects HR:partyRepublican` = 0"))
temp_df_2 <- 
  as.data.frame(matrix(c("CATE on Republicans", 
                         confint(temp_estimate_2)$confint[ , c("Estimate", "lwr", "upr")], 
                         "attack", "Ally,\nProtects HR"), 
                       nrow = 1))
temp_df_2 <- temp_df_2 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)


temp_estimate_3 <- glht(mod, linfct = c("`treatment_catNon-Ally, Protects HR` + `treatment_catNon-Ally, Protects HR:partyRepublican` = 0"))
temp_df_3 <- 
  as.data.frame(matrix(c("CATE on Republicans", 
                         confint(temp_estimate_3)$confint[ , c("Estimate", "lwr", "upr")], 
                         "attack", "Non-Ally,\nProtects HR"), 
                       nrow = 1))
temp_df_3 <- temp_df_3 %>% 
  rename(term = V1, estimate = V2, conf.low = V3, conf.high = V4, outcome = V5, treatment = V6)




CATE_rep_only <- bind_rows(temp_df_1, temp_df_2, temp_df_3) %>%
  mutate(term = as.factor(term),
         estimate = as.numeric(estimate),
         conf.high = as.numeric(conf.high),
         conf.low = as.numeric(conf.low),
         outcome = as.factor(outcome),
         treatment_cat = as.factor(treatment))

# Join together
CATE_total_df <- CATE_summary_df %>%
  mutate(conf.high = estimate + (1.96*std.error),
         conf.low  = estimate - (1.96*std.error))%>%
  bind_rows(CATE_rep_only)%>%
  mutate(term = factor(term,
                       levels = c("CATE on Democrats", 
                                  "CATE on Republicans",
                                  "Difference in CATEs")),
         treatment_cat = factor(treatment_cat,
                                levels = c("Ally,\nViolates HR",
                                           "Non-Ally,\nProtects HR",
                                           "Ally,\nProtects HR")))

###Plot Figure S1

p <- ggplot(CATE_total_df, aes(x=treatment_cat, y=estimate, color = term, shape = term)) + 
  geom_point(position = position_dodge(.5), size = 3) +
  #  scale_color_manual("", values = c("grey0", "grey40", "grey80")) +
  scale_color_manual("", values = c("blue", "red", "grey0")) +
  scale_shape_manual("", values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(.5))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  xlab("") + 
  ylab("Estimated Treatment Effect on\nSupport for Attack") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 12, family = "Times"),
        axis.text = element_text(color = "black", family = "Times", size = 11),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.key = element_rect(fill = "white"),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) 

p
# ggsave("CATE_party.tiff", width = 10, height = 5)
```

```{r}
#### Figure S2: Marginal Treatment Effect of Support for Attack by Nationalism ####
#Order of treatments
#Ally, Violates HR
#Non-Ally, Protects HR
#Ally, Protects HR

#Create interaction model between treatment and nationalism
mod_nat <- lm(attack ~ treatment_cat*nationalism_rc + male + age_cat + edu4 + inc + white, data = df_heter)

df_plot <- plot_slopes(mod_nat, 
                       variables = "treatment_cat", 
                       condition = "nationalism_rc",
                       conf_level = 0.95,
                       draw = FALSE) %>%
  #rename for plotting  
  mutate(contrast = ifelse(contrast == "Ally, Protects HR - Non-Ally, Violates HR", "Ally, Protects HR",
                           ifelse(contrast == "Ally, Violates HR - Non-Ally, Violates HR", "Ally, Violates HR", 
                                  ifelse(contrast == "Non-Ally, Protects HR - Non-Ally, Violates HR", "Non-Ally, Protects HR", contrast))),
         contrast = factor(contrast,
                           levels = c("Ally, Violates HR",
                                      "Non-Ally, Protects HR",
                                      "Ally, Protects HR")))

#Plot heterogeneous effects
p <- ggplot(df_plot, aes(y = estimate, x = nationalism_rc)) +
  geom_line() +
  facet_wrap(facets = vars(contrast)) +
  geom_ribbon(width = 0, aes(ymin = conf.low, ymax = conf.high), alpha = .3)+
  ylim(-30, 15) + 
  geom_hline(yintercept=0, linetype="dashed", color = "black")+
  labs(y = "Marginal Treatment Effect on\nSupport for Attack (95% CIs)",
       x = "Nationalism (1 = least nationalistic; 5 = most nationalistic)")+
  theme_bw() +
  theme(text = element_text(color = "black", size = 12, family = "Times"),
        axis.text = element_text(color = "black", family = "Times", size = 11))
p
# ggsave("het_effects_nat.tiff", width = 10, height = 5)

```


```{r}
#### Figure S3: Marginal Treatment Effect of Support for Attack by Patriotism ####

#Create interaction model between treatment and patriotism
mod_pat <- lm(attack ~ treatment_cat* patriotism_rc + male + age_cat + edu4 + inc + white, data = df_heter)

df_plot_pat <- plot_slopes(mod_pat, 
                           variables = "treatment_cat", 
                           condition = "patriotism_rc",
                           conf_level = 0.95,
                           draw = FALSE) %>%
  #rename for plotting  
  mutate(contrast = ifelse(contrast == "Ally, Protects HR - Non-Ally, Violates HR", "Ally, Protects HR",
                           ifelse(contrast == "Ally, Violates HR - Non-Ally, Violates HR", "Ally, Violates HR", 
                                  ifelse(contrast == "Non-Ally, Protects HR - Non-Ally, Violates HR", "Non-Ally, Protects HR", contrast))),
         contrast = factor(contrast,
                           levels = c("Ally, Violates HR",
                                      "Non-Ally, Protects HR",
                                      "Ally, Protects HR")))

#Plot heterogeneous effects
p <- ggplot(df_plot_pat, aes(y = estimate, x = patriotism_rc)) +
  geom_line() +
  facet_wrap(facets = vars(contrast)) +
  geom_ribbon(width = 0, aes(ymin = conf.low, ymax = conf.high), alpha = .3)+
  ylim(-30, 15) + 
  geom_hline(yintercept=0, linetype="dashed", color = "black")+
  labs(y = "Marginal Treatment Effect on\nSupport for Attack (95% CIs)",
       x = "Patriotism (1 = least patriotic; 4 = most patriotic)")+
  theme_bw() +
  theme(text = element_text(color = "black", size = 12, family = "Times"),
        axis.text = element_text(color = "black", family = "Times", size = 11))
p
# ggsave("het_effects_pat.tiff", width = 10, height = 5)
```

```{r}
#### Figure S4: Marginal Treatment Effect of Support for Attack by Cooperative Internationalism ####

#Create interaction model between treatment and coop_int
mod_coop <- lm(attack ~ treatment_cat*coop_int + male + age_cat + edu4 + inc + white, data = df_heter)

df_plot_coop <- plot_slopes(mod_coop, 
                            variables = "treatment_cat", 
                            condition = "coop_int",
                            conf_level = 0.95,
                            draw = FALSE) %>%
  #rename for plotting  
  mutate(contrast = ifelse(contrast == "Ally, Protects HR - Non-Ally, Violates HR", "Ally, Protects HR",
                           ifelse(contrast == "Ally, Violates HR - Non-Ally, Violates HR", "Ally, Violates HR", 
                                  ifelse(contrast == "Non-Ally, Protects HR - Non-Ally, Violates HR", "Non-Ally, Protects HR", contrast))),
         contrast = factor(contrast,
                           levels = c("Ally, Violates HR",
                                      "Non-Ally, Protects HR",
                                      "Ally, Protects HR")))

#Plot heterogeneous effects
p <- ggplot(df_plot_coop, aes(y = estimate, x = coop_int)) +
  geom_line() +
  facet_wrap(facets = vars(contrast)) +
  geom_ribbon(width = 0, aes(ymin = conf.low, ymax = conf.high), alpha = .3)+
  ylim(-30, 15) + 
  geom_hline(yintercept=0, linetype="dashed", color = "black")+
  labs(y = "Marginal Treatment Effect on\nSupport for Attack (95% CIs)",
       x = "Cooperative Internationalism (1 = least cooperative; 5 = most cooperative)")+
  theme_bw() +
  theme(text = element_text(color = "black", size = 12, family = "Times"),
        axis.text = element_text(color = "black", family = "Times", size = 11))
p
# ggsave("het_effects_coop.tiff", width = 10, height = 5)
```

```{r}
#### Figure S5: Mechanisms as Percentage of Total Effect ####

### (S5a) Replot of Tomz and Weeks (2020)
# the estimates are calculated from original dataset of Tomz and Weeks
prop_threat_tw <- 26.476151
threat_ci_low_tw <- 16.97482
threat_ci_high_tw <- 35.97748

prop_moral_tw <- 40.518353   
moral_ci_low_tw <- 28.09105
moral_ci_high_tw <- 52.94566

prop_success_tw <- 1.5941233  
success_ci_low_tw <- -.8727255   
success_ci_high_tw <- 4.060972

prop_cost_tw <- 1.6944775
cost_ci_low_tw <- -1.095952
cost_ci_high_tw <- 4.484907

prop_other_tw <- 29.716895
other_ci_low_tw <- 12.07185
other_ci_high_tw <- 47.36194

# aggregate the estimates in a dataframe and plot
fig4_tw_r <- data.frame(y = c(prop_threat_tw, prop_moral_tw, prop_success_tw, prop_cost_tw, prop_other_tw), 
                        x = as.factor(c("Threat", "Moral", "Success", "Cost", "Other")),
                        ci_low = c(threat_ci_low_tw, moral_ci_low_tw, success_ci_low_tw, cost_ci_low_tw, other_ci_low_tw),
                        ci_high = c(threat_ci_high_tw, moral_ci_high_tw, success_ci_high_tw, cost_ci_high_tw, other_ci_high_tw))
positions <- c("Threat", "Moral", "Success", "Cost", "Other")

p <- ggplot(fig4_tw_r, aes(x=x, y=y, ymin=ci_low, ymax=ci_high))+ 
  scale_x_discrete(limits = rev(positions)) +
  geom_pointrange()+
  coord_flip()+ 
  geom_text(aes(label=round(y, 0)), position=position_dodge(width=0.9), vjust=-0.75) +
  xlab('Mediator') +
  ylab("Effect of Human Rights on Support for War") + 
  theme_classic()
p
# ggsave("fig4_tw_r.pdf", width = 4.5, height = 3)

### (S5b) Replication Using PureSpectrum Survey

# use the product of regression coefficients method by Baron and Kenney (1986)
# first calculate the effect of mediators on DV
threat <- lm(threat ~ alliance*hmrts + male + age + edu4 + income, data = df_clean)
threat_on_DV <- as.numeric(threat$coef["hmrts"] + threat$coef["alliance:hmrts"])

moral <- lm(moral ~ alliance*hmrts + male + age + edu4 + income, data = df_clean)
moral_on_DV <- as.numeric(moral$coef["hmrts"] + moral$coef["alliance:hmrts"])

success <- lm(success ~ alliance*hmrts + male + age + edu4 + income, data = df_clean)
success_on_DV <- as.numeric(success$coef["hmrts"] + success$coef["alliance:hmrts"])

cost <- lm(cost ~ alliance*hmrts + male + age + edu4 + income, data = df_clean)
cost_on_DV <- as.numeric(cost$coef["hmrts"] + cost$coef["alliance:hmrts"])

# second calculate the effect of treatment on mediators
dv_on_med <- lm(attack ~ threat + moral + success + cost + alliance*hmrts + 
                  male + age + edu4 + income, data = df_clean)

dv_on_threat <- as.numeric(dv_on_med$coef["threat"])
dv_on_moral <- as.numeric(dv_on_med$coef["moral"])
dv_on_success <- as.numeric(dv_on_med$coef["success"])
dv_on_cost <- as.numeric(dv_on_med$coef["cost"])

# third calculate the total treatment effect
total <- lm(attack ~ alliance*hmrts + male + age + edu4 + income, data = df_clean)
total_dv <- as.numeric(total$coef["hmrts"] + total$coef["alliance:hmrts"])

# finally calculate the proportion explained
prop_threat <- (threat_on_DV*dv_on_threat)/total_dv
prop_moral <- (moral_on_DV*dv_on_moral)/total_dv
prop_success <- (success_on_DV*dv_on_success)/total_dv
prop_cost <- (cost_on_DV*dv_on_cost)/total_dv
prop_other <- 1 - prop_threat - prop_moral - prop_success - prop_cost

# bootstrap CIs on proportion explained
set.seed(1234)

# threat
threat_boot <- function(dat, i){
  df <- dat[i,]
  
  threat.fit <- lm(threat ~ alliance*hmrts + male + age + edu4 + income, data = df)
  threat_on_DV <- as.numeric(threat.fit$coef["hmrts"] + threat.fit$coef["alliance:hmrts"])
  
  dv_on_med <- lm(attack ~ threat + moral + success + cost + alliance*hmrts + 
                    male + age + edu4 + income, data = df)
  dv_on_threat <- as.numeric(dv_on_med$coef["threat"])
  
  total <- lm(attack ~ alliance*hmrts + male + age + edu4 + income, data = df)
  total_dv <- as.numeric(total$coef["hmrts"] + total$coef["alliance:hmrts"])
  
  (threat_on_DV*dv_on_threat)/total_dv
}

boot_threat <- boot(df_clean, threat_boot, 1000, stype = "i")
threat_ci_low <- prop_threat - 1.96*sd(boot_threat$t)
threat_ci_high <- prop_threat + 1.96*sd(boot_threat$t)

# moral
moral_boot <- function(dat, i){
  df <- dat[i,]
  
  moral.fit <- lm(moral ~ alliance*hmrts + male + age + edu4 + income, data = df)
  moral_on_DV <- as.numeric(moral.fit$coef["hmrts"] + moral.fit$coef["alliance:hmrts"])
  
  dv_on_med <- lm(attack ~ threat + moral + success + cost + alliance*hmrts + 
                    male + age + edu4 + income, data = df)
  dv_on_moral <- as.numeric(dv_on_med$coef["moral"])
  
  total <- lm(attack ~ alliance*hmrts + male + age + edu4 + income, data = df)
  total_dv <- as.numeric(total$coef["hmrts"] + total$coef["alliance:hmrts"])
  
  (moral_on_DV*dv_on_moral)/total_dv
}
boot_moral <- boot(df_clean, moral_boot, 1000, stype = "i")
moral_ci_low <- prop_moral - 1.96*sd(boot_moral$t)
moral_ci_high <- prop_moral + 1.96*sd(boot_moral$t)

# success
success_boot <- function(dat, i){
  df <- dat[i,]
  
  success.fit <- lm(success ~ alliance*hmrts + male + age + edu4 + income, data = df)
  success_on_DV <- as.numeric(success.fit$coef["hmrts"] + success.fit$coef["alliance:hmrts"])
  
  dv_on_med <- lm(attack ~ threat + moral + success + cost + alliance*hmrts + 
                    male + age + edu4 + income, data = df)
  dv_on_success <- as.numeric(dv_on_med$coef["success"])
  
  total <- lm(attack ~ alliance*hmrts + male + age + edu4 + income, data = df)
  total_dv <- as.numeric(total$coef["hmrts"] + total$coef["alliance:hmrts"])
  
  (success_on_DV*dv_on_success)/total_dv
}
boot_success <- boot(df_clean, success_boot, 1000, stype = "i")
success_ci_low <- prop_success - 1.96*sd(boot_success$t)
success_ci_high <- prop_success + 1.96*sd(boot_success$t)

# cost
cost_boot <- function(dat, i){
  df <- dat[i,]
  
  cost.fit <- lm(cost ~ alliance*hmrts + male + age + edu4 + income, data = df)
  cost_on_DV <- as.numeric(cost.fit$coef["hmrts"] + cost.fit$coef["alliance:hmrts"])
  
  dv_on_med <- lm(attack ~ threat + moral + success + cost + alliance*hmrts + 
                    male + age + edu4 + income, data = df)
  dv_on_cost <- as.numeric(dv_on_med$coef["cost"])
  
  total <- lm(attack ~ alliance*hmrts + male + age + edu4 + income, data = df)
  total_dv <- as.numeric(total$coef["hmrts"] + total$coef["alliance:hmrts"])
  
  (cost_on_DV*dv_on_cost)/total_dv
}
boot_cost <- boot(df_clean, cost_boot, 1000, stype = "i")
cost_ci_low <- prop_cost - 1.96*sd(boot_cost$t)
cost_ci_high <- prop_cost + 1.96*sd(boot_cost$t)

# other
other_boot <- function(dat, i){
  df <- dat[i,]
  threat.fit <- lm(threat ~ alliance*hmrts + male + age + edu4 + income, data = df)
  threat_on_DV <- as.numeric(threat.fit$coef["hmrts"] + threat.fit$coef["alliance:hmrts"])
  
  moral.fit <- lm(moral ~ alliance*hmrts + male + age + edu4 + income, data = df)
  moral_on_DV <- as.numeric(moral.fit$coef["hmrts"] + moral.fit$coef["alliance:hmrts"])
  
  success.fit <- lm(success ~ alliance*hmrts + male + age + edu4 + income, data = df)
  success_on_DV <- as.numeric(success.fit$coef["hmrts"] + success.fit$coef["alliance:hmrts"])
  
  cost.fit <- lm(cost ~ alliance*hmrts + male + age + edu4 + income, data = df)
  cost_on_DV <- as.numeric(cost.fit$coef["hmrts"] + cost.fit$coef["alliance:hmrts"])
  
  dv_on_med <- lm(attack ~ threat + moral + success + cost + alliance*hmrts + 
                    male + age + edu4 + income, data = df)
  dv_on_threat <- as.numeric(dv_on_med$coef["threat"])
  dv_on_moral <- as.numeric(dv_on_med$coef["moral"])
  dv_on_success <- as.numeric(dv_on_med$coef["success"])
  dv_on_cost <- as.numeric(dv_on_med$coef["cost"])
  
  total <- lm(attack ~ alliance*hmrts + male + age + edu4 + income, data = df)
  total_dv <- as.numeric(total$coef["hmrts"] + total$coef["alliance:hmrts"])
  
  prop_threat <- (threat_on_DV*dv_on_threat)/total_dv
  prop_moral <- (moral_on_DV*dv_on_moral)/total_dv
  prop_success <- (success_on_DV*dv_on_success)/total_dv
  prop_cost <- (cost_on_DV*dv_on_cost)/total_dv
  1 - prop_threat - prop_moral - prop_success - prop_cost
}
boot_other <- boot(df_clean, other_boot, 1000, stype = "i")
other_ci_low <- prop_other - 1.96*sd(boot_other$t)
other_ci_high <- prop_other + 1.96*sd(boot_other$t)


# plot proportion estimated
prop_est_hmrts <- data.frame(y = c(prop_threat, prop_moral, prop_success, prop_cost, prop_other), 
                             x = as.factor(c("Threat", "Moral", "Success", "Cost", "Other")),
                             ci_low = c(threat_ci_low, moral_ci_low, success_ci_low, cost_ci_low, other_ci_low),
                             ci_high = c(threat_ci_high, moral_ci_high, success_ci_high, cost_ci_high, other_ci_high))
positions <- c("Threat", "Moral", "Success", "Cost", "Other")
p1 <- ggplot(prop_est_hmrts, aes(x=x, y=y*100, ymin=ci_low*100, ymax=ci_high*100))+ 
  scale_x_discrete(limits = rev(positions)) +
  geom_pointrange()+
  coord_flip()+ 
  geom_text(aes(label=round(y, 2)*100), position=position_dodge(width=0.9), vjust=-0.75) +
  xlab('Mediator') +
  ylab("Effect of Human Rights on Support for War") + 
  theme_classic()

p1
# ggsave("prop_est_hmrts.pdf", width = 4.5, height = 3)

```

```{r}
#### Table S4: Causal Mediation Analysis (Potential Outcomes Framework) ####

## calculate the ACMEs of the Tomz and Weeks data
# load the Tomz and Weeks replication data
df_tw <- read_dta("2012-10-01-Main-prepped.dta")
#recode TW hrts variable so 1 = violates human rights
df_tw <- df_tw %>% mutate(hrts = ifelse(hrts == 0, 1, 0)) 

# 100 repetitions of n=100 simulations for efficiency
sims.n<-100
reps.n<-100

# calculate ACME for each mediator
set.seed(12345)

## THREATS ##
out.t <-lm(strike ~ democ*hrts + threat + 
             h1 + i1 + p1 + e1 + r1 + male + white + age + ed4, data=df_tw) 
med.t <- lm(threat ~ democ*hrts + h1 + i1 + p1 + e1 + r1 + male + white + age + ed4, data=df_tw)
med.out.t <- mediate(med.t, out.t, treat="hrts", mediator="threat", sims=100)

prop.med.t<-c()
acme.med.t<-c()
for(i in 1:reps.n){
  cma.t <- mediate(med.t, out.t, treat="hrts", mediator="threat", sims=sims.n)
  prop.med.t[i]<-summary(cma.t)$n.avg
  acme.med.t[i]<-summary(cma.t)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## MORAL ##
out.m <-lm(strike ~ democ*hrts + moral + 
             h1 + i1 + p1 + e1 + r1 + male + white + age + ed4, data=df_tw) 
med.m <- lm(moral ~ democ*hrts + h1 + i1 + p1 + e1 + r1 + male + white + age + ed4, data=df_tw)
med.out.m <- mediate(med.m, out.m, treat="hrts", mediator="moral", sims=100)

prop.med.m<-c()
acme.med.m<-c()
i <- 1
for(i in 1:reps.n){
  cma.m <- mediate(med.m, out.m, treat="hrts", mediator="moral", sims=sims.n)
  prop.med.m[i]<-summary(cma.m)$n.avg
  acme.med.m[i]<-summary(cma.m)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## SUCCESS ##
out.s <-lm(strike ~ democ*hrts + success + 
             h1 + i1 + p1 + e1 + r1 + male + white + age + ed4, data=df_tw) 
med.s <- lm(success ~ democ*hrts + h1 + i1 + p1 + e1 + r1 + male + white + age + ed4, data=df_tw)
med.out.s <- mediate(med.s, out.s, treat="hrts", mediator="success", sims=100)

prop.med.s<-c()
acme.med.s<-c()
i <- 1
for(i in 1:reps.n){
  cma.s <- mediate(med.s, out.s, treat="hrts", mediator="success", sims=sims.n)
  prop.med.s[i]<-summary(cma.s)$n.avg
  acme.med.s[i]<-summary(cma.s)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## COST ##
out.c <-lm(strike ~ democ*hrts + cost + 
             h1 + i1 + p1 + e1 + r1 + male + white + age + ed4, data=df_tw) 
med.c <- lm(cost ~ democ*hrts + h1 + i1 + p1 + e1 + r1 + male + white + age + ed4, data=df_tw)
med.out.c <- mediate(med.c, out.c, treat="hrts", mediator="cost", sims=100)

prop.med.c<-c()
acme.med.c<-c()
i <- 1
for(i in 1:reps.n){
  cma.c <- mediate(med.c, out.c, treat="hrts", mediator="cost", sims=sims.n)
  prop.med.c[i]<-summary(cma.c)$n.avg
  acme.med.c[i]<-summary(cma.c)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## calculate the ACMEs of the PureSpectrum data using Human Rights Violation as treatment

## THREATS ##
out.t.rep <-lm(attack ~ threat + hmrts*alliance + male + edu4 + age + income, data=df_clean) 
med.t.rep <- lm(threat ~ hmrts*alliance + male + edu4 + age + income, data=df_clean)
med.out.t.rep <- mediate(med.t.rep, out.t.rep, treat="hmrts", mediator="threat", sims=100)

prop.med.t.rep<-c()
acme.med.t.rep<-c()
i <- 1
for(i in 1:reps.n){
  cma.t.rep <- mediate(med.t.rep, out.t.rep, treat="hmrts", mediator="threat", sims=sims.n)
  prop.med.t.rep[i]<-summary(cma.t.rep)$n.avg
  acme.med.t.rep[i]<-summary(cma.t.rep)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## MORAL ##
out.m.rep <-lm(attack ~ moral + hmrts*alliance + male + edu4 + age + income, data=df_clean) 
med.m.rep <- lm(moral ~ hmrts*alliance + male + edu4 + age + income, data=df_clean)
med.out.m.rep <- mediate(med.m.rep, out.m.rep, treat="hmrts", mediator="moral", sims=100)

prop.med.m.rep<-c()
acme.med.m.rep<-c()
i <- 1
for(i in 1:reps.n){
  cma.m.rep <- mediate(med.m.rep, out.m.rep, treat="hmrts", mediator="moral", sims=sims.n)
  prop.med.m.rep[i]<-summary(cma.m.rep)$n.avg
  acme.med.m.rep[i]<-summary(cma.m.rep)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## SUCCESS ##
out.s.rep <-lm(attack ~ success + hmrts*alliance + male + edu4 + age + income, data=df_clean) 
med.s.rep <- lm(success ~ hmrts*alliance + male + edu4 + age + income, data=df_clean)
med.out.s.rep <- mediate(med.s.rep, out.s.rep, treat="hmrts", mediator="success", sims=100)

prop.med.s.rep<-c()
acme.med.s.rep<-c()
i <- 1
for(i in 1:reps.n){
  cma.s.rep <- mediate(med.s.rep, out.s.rep, treat="hmrts", mediator="success", sims=sims.n)
  prop.med.s.rep[i]<-summary(cma.s.rep)$n.avg
  acme.med.s.rep[i]<-summary(cma.s.rep)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## COST ##
out.c.rep <-lm(attack ~ cost + hmrts*alliance + male + edu4 + age + income, data=df_clean) 
med.c.rep <- lm(cost ~ hmrts*alliance + male + edu4 + age + income, data=df_clean)
med.out.c.rep <- mediate(med.c.rep, out.c.rep, treat="hmrts", mediator="cost", sims=100)

prop.med.c.rep<-c()
acme.med.c.rep<-c()
i <- 1
for(i in 1:reps.n){
  cma.c.rep <- mediate(med.c.rep, out.c.rep, treat="hmrts", mediator="cost", sims=sims.n)
  prop.med.c.rep[i]<-summary(cma.c.rep)$n.avg
  acme.med.c.rep[i]<-summary(cma.c.rep)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## calculate the ACMEs of the PureSpectrum data using U.S. Alliance as treatment
## THREATS ##
out.t.alliance <-lm(attack ~ threat + hmrts*alliance + male + edu4 + age + income, data=df_clean) 
med.t.alliance <- lm(threat ~ hmrts*alliance + male + edu4 + age + income, data=df_clean)
med.out.t.alliance <- mediate(med.t.alliance, out.t.alliance, treat="alliance", mediator="threat", sims=100)

prop.med.t.alliance<-c()
acme.med.t.alliance<-c()
i <- 1
for(i in 1:reps.n){
  cma.t.alliance <- mediate(med.t.alliance, out.t.alliance, treat="alliance", mediator="threat", sims=sims.n)
  prop.med.t.alliance[i]<-summary(cma.t.alliance)$n.avg
  acme.med.t.alliance[i]<-summary(cma.t.alliance)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## MORAL ##
out.m.alliance <-lm(attack ~ moral + hmrts*alliance + male + edu4 + age + income, data=df_clean) 
med.m.alliance <- lm(moral ~ hmrts*alliance + male + edu4 + age + income, data=df_clean)
med.out.m.alliance <- mediate(med.m.alliance, out.m.alliance, treat="alliance", mediator="moral", sims=100)

prop.med.m.alliance<-c()
acme.med.m.alliance<-c()
i <- 1
for(i in 1:reps.n){
  cma.m.alliance <- mediate(med.m.alliance, out.m.alliance, treat="alliance", mediator="moral", sims=sims.n)
  prop.med.m.alliance[i]<-summary(cma.m.alliance)$n.avg
  acme.med.m.alliance[i]<-summary(cma.m.alliance)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## SUCCESS ##
out.s.alliance <-lm(attack ~ success + hmrts*alliance + male + edu4 + age + income, data=df_clean) 
med.s.alliance <- lm(success ~ hmrts*alliance + male + edu4 + age + income, data=df_clean)
med.out.s.alliance <- mediate(med.s.alliance, out.s.alliance, treat="alliance", mediator="success", sims=100)

prop.med.s.alliance<-c()
acme.med.s.alliance<-c()
i <- 1
for(i in 1:reps.n){
  cma.s.alliance <- mediate(med.s.alliance, out.s.alliance, treat="alliance", mediator="success", sims=sims.n)
  prop.med.s.alliance[i]<-summary(cma.s.alliance)$n.avg
  acme.med.s.alliance[i]<-summary(cma.s.alliance)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## COST ##
out.c.alliance <-lm(attack ~ cost + hmrts*alliance + male + edu4 + age + income, data=df_clean) 
med.c.alliance <- lm(cost ~ hmrts*alliance + male + edu4 + age + income, data=df_clean)
med.out.c.alliance <- mediate(med.c.alliance, out.c.alliance, treat="alliance", mediator="cost", sims=100)

prop.med.c.alliance<-c()
acme.med.c.alliance<-c()
i <- 1
for(i in 1:reps.n){
  cma.c.alliance <- mediate(med.c.alliance, out.c.alliance, treat="alliance", mediator="cost", sims=sims.n)
  prop.med.c.alliance[i]<-summary(cma.c.alliance)$n.avg
  acme.med.c.alliance[i]<-summary(cma.c.alliance)$d.avg
  if(i %% 10 == 0)	{ cat(i,"\n") }
}

## compile the table
tab1 <- data.frame(matrix(ncol = 3, nrow = 4))
row.names(tab1) <- c("Threat", "Moral", "Success", "Cost")

tab1[,1] <- c(round(mean(acme.med.t),2), round(mean(acme.med.m),2), 
              round(mean(acme.med.s),2), round(mean(acme.med.c),2))
tab1[,2] <- c(round(mean(acme.med.t.rep),2), round(mean(acme.med.m.rep),2), 
              round(mean(acme.med.s.rep),2), round(mean(acme.med.c.rep),2))
tab1[,3] <- c(round(mean(acme.med.t.alliance),2), round(mean(acme.med.m.alliance),2), 
              round(mean(acme.med.s.alliance),2), round(mean(acme.med.c.alliance),2))
tab1

```

```{r}
#### Figure S6: Mediation Sensitivity with Respect to Different Sensitivity Parameter \rho ####

# conducting sensitivity analysis for the mediation analysis in Table S4

## THREATS ##
sens.out.t <- medsens(med.out.t, rho.by = 0.1, effect.type = "indirect", sims = 100)

## MORAL ##
sens.out.m <- medsens(med.out.m, rho.by = 0.1, effect.type = "indirect", sims = 100)

## SUCCESS ##
sens.out.s <- medsens(med.out.s, rho.by = 0.1, effect.type = "indirect", sims = 100)

## COST ##
sens.out.c <- medsens(med.out.c, rho.by = 0.1, effect.type = "indirect", sims = 100)

## PureSpectrum data using Human Rights Violation as treatment

## THREATS ##
sens.out.t.rep <- medsens(med.out.t.rep, rho.by = 0.1, effect.type = "indirect", sims = 100)

## MORAL ##
sens.out.m.rep <- medsens(med.out.m.rep, rho.by = 0.1, effect.type = "indirect", sims = 100)

## SUCCESS ##
sens.out.s.rep <- medsens(med.out.s.rep, rho.by = 0.1, effect.type = "indirect", sims = 100)

## COST ##
sens.out.c.rep <- medsens(med.out.c.rep, rho.by = 0.1, effect.type = "indirect", sims = 100)

## PureSpectrum data using U.S. Alliance as treatment

## THREATS ##
sens.out.t.alliance <- medsens(med.out.t.alliance, rho.by = 0.1, effect.type = "indirect", sims = 100)

## MORAL ##
sens.out.m.alliance <- medsens(med.out.m.alliance, rho.by = 0.1, effect.type = "indirect", sims = 100)

## SUCCESS ##
sens.out.s.alliance <- medsens(med.out.s.alliance, rho.by = 0.1, effect.type = "indirect", sims = 100)

## COST ##
sens.out.c.alliance <- medsens(med.out.c.alliance, rho.by = 0.1, effect.type = "indirect", sims = 100)

# aggregate the sensitivity analysis plots on the same page
# create three long plots, grouped by survey and treatment

# TW data, Tr = Human Rights Violation
# pdf(file= "sensitivity1.pdf",width = 4, height = 12)
# par(mfrow = c(4, 1))
# TW, threat
plot(sens.out.t, sens.par = "rho", main = "(a1) Threat",
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")

# TW, moral
plot(sens.out.m, sens.par = "rho", main = "(a2) Moral", 
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")

# TW, success
plot(sens.out.s, sens.par = "rho", main = "(a3) Success",
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")

# TW, cost
plot(sens.out.c, sens.par = "rho", main = "(a4) Cost", 
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")
# dev.off()

# PureSpectrum data, Tr = Human Rights Violation
# pdf(file= "sensitivity2.pdf",width = 4, height = 12)
# par(mfrow = c(4, 1))
#PS, tr=hmrts, threat
plot(sens.out.t.rep, sens.par = "rho", main = "(b1) Threat", 
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")

#PS, tr=hmrts, moral
plot(sens.out.m.rep, sens.par = "rho", main = "(b2) Moral", 
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")

#PS, tr=hmrts, success
plot(sens.out.s.rep, sens.par = "rho", main = "(b3) Success", 
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")

#PS, tr=hmrts, cost
plot(sens.out.c.rep, sens.par = "rho", main = "(b4) Cost", 
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")
# dev.off()

# PureSpectrum data, Tr = U.S. Alliance
# pdf(file= "sensitivity3.pdf",width = 4, height = 12)
# par(mfrow = c(4, 1))
#PS, tr=alliance, threat
plot(sens.out.t.alliance, sens.par = "rho", main = "(c1) Threat", 
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")

#PS, tr=alliance, moral
plot(sens.out.m.alliance, sens.par = "rho", main = "(c2) Moral", 
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")

#PS, tr=alliance, success
plot(sens.out.s.alliance, sens.par = "rho", main = "(c3) Success", 
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")

#PS, tr=alliance, cost
plot(sens.out.c.alliance, sens.par = "rho", main = "(c4) Cost", 
     xlab = expression(paste("Sensitivity Parameter ", rho)), ylab = "Average Mediation Effect")
# dev.off()

```

```{r, out.extra='angle=90', results='asis'}
#### Table S5: Regression Estimates of Support for War (5-Point Likert Scale, PureSpectrum Survey) ####

# run the regression models
m1_cont <- lm(attack_cont ~ hmrts + alliance, data = df_clean)
m2_cont <- lm(attack_cont ~ alliance * hmrts, data = df_clean)
m3_cont <- lm(attack_cont ~ alliance*hmrts + 
                male + age_cat + edu4 + inc_10k, data = df_clean)

texreg(l = list(m1_cont, m2_cont, m3_cont),
       reorder.coef= c(2, 3, 4, 5, 6, 7, 8, 1),
       custom.coef.names = c("(Intercept)", "Violating Human Rights", "U.S. Military Alliance",
                             "Violating Human Rights $\\times$ U.S. Military Alliance",
                             "Male", "Age", "Education", "Income"),
       stars = c(0.01, 0.05, 0.1),
       digits = 2,
       caption = "Regression Estimates of Support for War (5-point Likert Scale, PureSpectrum Survey)",
       caption.above = T,
       include.ci = F,
       include.rmse =F,
       include.rsq = F,
       include.adjrs = F,
       label = "",
       custom.note = "",
       fontsize = "small") %>%
  gsub(".begin.center.", "\\\\centering", .) %>%
  gsub(".end.center.", "", .)

```

```{r}
#### Figure S7: Impact of Treatments on Support for War (5-Point Likert Scale, PureSpectrum Survey, 95% Confidence Intervals) ####

### Plot (S7a) Support for Attack
# calculate the mean support for attack in each treatment group
df_ate_cont <- df_clean %>% 
  group_by(exp_4, alliance, hmrts) %>%
  summarise(ate = mean(attack_cont, na.rm = TRUE),
            n = n(),
            sd = sd(attack, na.rm = TRUE),
            se = sd(attack, na.rm = TRUE) / sqrt(n)) %>%
  mutate(ci_low = ate - 1.96*se,
         ci_high = ate + 1.96*se)

p <- ggplot(df_ate_cont, aes(x = factor(exp_4, level=c('1', '3', '2', '4')), y = ate,
                        color=factor(exp_4))) + 
  theme_classic() +
  geom_point()+
  geom_errorbar(aes(ymin=ci_low, ymax=ci_high), width=.1,
                position=position_dodge(0.05)) +
  scale_x_discrete(labels= c('Non-Ally,Violates HR', 'Ally,Violates HR',
                             'Non-Ally,Protects HR', "Ally,Protects HR")) +
  scale_color_manual(values=c('#999999', 'black', 'black', 'black')) +
  theme(legend.position = "none") +
  labs(x = "", y = "Mean Support for Attack \n (% point)", size = 12) +
  # theme(axis.title.y = element_text(size = 9)) +
  geom_text(aes(label=round(ate, 1)), position=position_dodge(width=0.9), 
            vjust=.5, hjust = -.35) +
  theme(axis.text.x = element_text(angle = 20, hjust = 0.5, vjust = 0.5, size = 12),
        axis.title.y = element_text(size=12))
p
# ggsave("ate_cont.pdf", width = 6, height = 4)

### Plot (S7b) Average Treatment Effect on Support for Attack
# calculate the difference in support for attack between 2-4 against baseline condition
est <- rep(NA, 4)
ci_low <- rep(NA, 4)
ci_high <- rep(NA, 4)
se <- rep(NA, 4)

for(i in 2:4){
  test <- t.test(df_clean$attack_cont[df_clean$exp_4==i], 
                 df_clean$attack_cont[df_clean$exp_4==1])
  est[i] <- test[["estimate"]][["mean of x"]] - test[["estimate"]][["mean of y"]]
  ci_low[i] <- test[["conf.int"]][1]
  ci_high[i] <- test[["conf.int"]][2]
  se[i] <- test[["stderr"]]
}

df_ate_cont_diff <- data.frame(exp_4 = df_ate_cont$exp_4, est, ci_low, ci_high, se)
df_ate_cont_diff <- df_ate_cont_diff[-1, ]

p1 <- ggplot(df_ate_cont_diff, aes(x = factor(exp_4, level=c('3', '2', '4')), y = est)) + 
  theme_classic() +
  geom_point()+
  geom_errorbar(aes(ymin=ci_low, ymax=ci_high), width=.1,
                position=position_dodge(0.05)) +
  scale_x_discrete(labels= c('Ally,Violates HR', 'Non-Ally,Protects HR', 
                             'Ally,Protects HR')) +
  labs(x = "", y = "Average Treatment Effect \n (% point)") +
  # theme(axis.title.y = element_text(size = 9)) +
  geom_text(aes(label=round(est, 1)), position=position_dodge(width=0.9), 
            vjust=.5, hjust = -.35) +
  geom_hline(yintercept = 0, linetype="dotted") +
  theme(axis.text.x = element_text(angle = 20, hjust = 0.5, vjust = 0.5, size = 12),
        axis.title.y = element_text(size=12))
p1
# ggsave("ate_cont-diff.pdf", width = 6, height = 4)

```

```{r, out.extra='angle=90', results='asis'}
#### Figure S8: The Estimated Effect of Alliance on Public Support for War (PureSpectrum Survey) ####

# reproduce the estimates
f <- attack ~ hmrts*alliance
fit <- lm(f, data = df_clean)

# the smallest substantively interesting effect
pure_high <- 4.3
pure_low <- 8.6

# a data frame setting the values of the "other" variables
X_c <- data.frame(
  alliance = 0,  # low value
  hmrts = 1 # low value
)

quantile(df_clean$alliance, c(.1, .9))
# 10th percentile of alliance = 0
# 90th percentile of alliance = 1

# compute the comparison for alliance
alliance_comp <- comparisons(fit,
                             newdata = X_c, 
                             variables = list("alliance" = c(0, 1)), # low to high value
                             conf_level = 0.90)

# bind the comparisons together and plot
comp <- bind_rows(alliance_comp, alliance_comp)
comp <- comp %>% dplyr::mutate(ID = row_number(),
                               study = c("PureSpectrum Survey", 
                                         "PureSpectrum Survey"),
                               threshold = c("A Stringent Threshold", "A Baseline Threshold"),
                               m = c(pure_high, pure_low),
                               label = c("m = 4.3", "m = 8.6"))

# plot the figure
set.seed(4755427)

# Create containers to store estimates and confidence intervals.
est.alliance <- matrix(NA, nrow = 1, ncol = 3)

# Plot our estimate
est.alliance[1, ] <- c(alliance_comp$conf.low, alliance_comp$estimate, alliance_comp$conf.high)
rownames(est.alliance) <- "The Estimate"

# pdf("equivalence-test-v2.pdf", height = 3, width = 9, family = "serif")
par(mfrow = c(1, 1), oma = c(.5, 9, .5, 9), mar = c(1,1,1,1))
eplot(NULL, xlim = c(-10, 10), ylim = c(-1, 0), 
      xat = c(-8.6, -4.3, -0.5, 4.3, 8.6),
      anny = FALSE,
      xlab = "",
      xlabpos = 2.5)
abline(v = 8.6, xpd = FALSE, lty = 2, col = "blue")
abline(v = -8.6, xpd = FALSE, lty = 2, col = "blue")
abline(v = -4.3, xpd = FALSE, lty = 5, col = "red")
abline(v = 4.3, xpd = FALSE, lty = 5, col = "red")
legend(par('usr')[2], par('usr')[4], bty='n', xpd=NA,
       c("Stringent Threshold", "Baseline Threshold"), lty=c(5,2), col = c("red", "blue"))

for (i in 1:nrow(est.alliance)) {
  est <- est.alliance
  lines(c(est[i, 1], est[i, 3]), c(-i/(nrow(est) + 1), -i/(nrow(est) + 1)), lwd = 2) 
  points(est[i, 2], -i/(nrow(est) + 1), pch = 19, cex = .7)
  text(est[i, 2], -i/(nrow(est) + 1), rownames(est)[i], pos = 3, cex = .8)
}
# dev.off()

```

```{r, out.extra='angle=90', results='asis'}
#### Table S6: Regression Estimates of Support for War (Binary Dependent Variable, Attentive Sample in PureSpectrum Survey) ####

# filter the attentive sample by checking if the respondents answered questions MVC_1 through MVC_3 correctly
df_atten <- df_clean %>% 
  filter(MVC_1 == 2 & MVC_2 == 1 & MVC_3 == 3)

# run the regression models
m1_atten <- lm(attack ~ hmrts+alliance, data = df_atten)
m2_atten <- lm(attack ~ alliance * hmrts, data = df_atten)
m3_atten <- lm(attack ~ alliance*hmrts + 
           male + age_cat + edu4 + inc_10k, data = df_atten)

texreg(l = list(m1_atten, m2_atten, m3_atten),
       reorder.coef= c(2, 3, 4, 5, 6, 7, 8, 1),
       custom.coef.names = c("(Intercept)", "Violating Human Rights", "U.S. Military Alliance",
                             "Violating Human Rights $\\times$ U.S. Military Alliance",
                             "Male", "Age", "Education", "Income"),
       stars = c(0.01, 0.05, 0.1),
       digits = 2,
       caption = "Regression Estimates of Support for War (Binary Dependent Variable, Attentive Sample in PureSpectrum Survey)",
       caption.above = T,
       include.ci = F,
       include.rmse =F,
       include.rsq = F,
       include.adjrs = F,
       custom.note = "",
       label = "tab:bin-dv",
       fontsize = "small") %>%
  gsub(".begin.center.", "\\\\centering", .) %>%
  gsub(".end.center.", "", .)


```

```{r, out.extra='angle=90', results='asis'}
#### Table S7: Regression Estimates of Support for War (5-Point Likert Scale, Attentive Sample in PureSpectrum Survey) ####

# run the regression models
m1_cont_atten <- lm(attack_cont ~ hmrts+alliance, data = df_atten)
m2_cont_atten <- lm(attack_cont ~ alliance * hmrts, data = df_atten)
m3_cont_atten <- lm(attack_cont ~ alliance*hmrts + 
                male + age_cat + edu4 + inc_10k, data = df_atten)

texreg(l = list(m1_cont_atten, m2_cont_atten, m3_cont_atten),
       reorder.coef= c(2, 3, 4, 5, 6, 7, 8, 1),
       custom.coef.names = c("(Intercept)", "Violating Human Rights", "U.S. Military Alliance",
                             "Violating Human Rights $\\times$ U.S. Military Alliance",
                             "Male", "Age", "Education", "Income"),
       stars = c(0.01, 0.05, 0.1),
       digits = 2,
       caption = "Regression Estimates of Support for War (Continuous Dependent Variable,
       Attentive Sample in PureSpectrum Survey)",
       caption.above = T,
       include.ci = F,
       include.rmse =F,
       include.rsq = F,
       include.adjrs = F,
       label = "",
       fontsize = "small") %>%
  gsub(".begin.center.", "\\\\centering", .) %>%
  gsub(".end.center.", "", .)

```

```{r}
#### Figure S9: Distribution of Pre-Treatment Variables (PureSpectrum Survey) ####

## recode the variables to be binary
var_balance <- df_clean %>% glimpse()

## female 
var_balance$female <- NA
var_balance$female [var_balance$sex==2]<-1 ## female
var_balance$female [var_balance$sex==1]<-0 ## male 

## white 
var_balance$white <- NA 
var_balance$white [var_balance$race==1]<-1 ## white
var_balance$white [var_balance$race>=2 & var_balance$race<=6]<-0 ## non-white  
table(var_balance$white)
## college degree 
var_balance$col<- NA
var_balance$col[var_balance$educ<=3]<-0 ## lower than 4-year college
var_balance$col[var_balance$educ>=4]<-1 ## college or above

## republican 
var_balance$republic<- NA
var_balance$republic [var_balance$pid_1==1] <- 1 ## republican 
var_balance$republic [var_balance$pid_1>=2]<-0 ## non-republican

## age
median(var_balance$age)
# median age 40-49
var_balance$midage <- NA
var_balance$midage[var_balance$age >= 3] <- 1 # age of 40 and above
var_balance$midage[var_balance$age <= 2] <- 0 # age of 39 and below

## ideology
var_balance$dem_eval <- NA
var_balance$dem_eval[var_balance$US_dem_eval  >= 6] <- 1 # democratic evaluation of 6 and above
var_balance$dem_eval[var_balance$US_dem_eval <= 5] <- 0 # democratic evaluation of 5 and below

## create the plot
# pdf(file= "var-balance.pdf",width = 7, height = 7)
# par(mfrow = c(2,3))
# par(mar = c(3.2,3.2,3.2,3.2))
# Gender
plot(density(var_balance[var_balance$exp_4==1,]$female, na.rm=T, bw=0.1), xlab="Gender (0 = Male, 1 = Female)", ylim=c(0,4), xlim=c(-0.4,1.4), 
     main="Gender \n(0 = Male, 1 = Female)", lty=1, cex=1.2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==2,]$female, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==3,]$female, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==4,]$female, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE, family="Times")
legend("topleft", legend=c("Not Ally, Violates HR", "Not Ally, Protects HR", "Ally, Violates HR", "Ally, Protects HR"), lty=c(1,2,3,4), cex=1, bty="n")


# Race
plot(density(var_balance[var_balance$exp_4==1,]$white, na.rm=T, bw=0.1), xlab="White (0 = No, 1 = Yes)", ylim=c(0,4), xlim=c(-0.4,1.4), 
     main="Race \n(0 = Non White, 1 = White)", lty=1, cex=1.2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==2,]$white, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==3,]$white, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==4,]$white, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE, family="Times")
legend("topleft", legend=c("Not Ally, Violates HR", "Not Ally, Protects HR", "Ally, Violates HR", "Ally, Protects HR"), lty=c(1,2,3,4), cex=1, bty="n")

# Education
plot(density(var_balance[var_balance$exp_4==1,]$col, na.rm=T, bw=0.1), xlab="College degree (0 = No, 1 = Yes)", ylim=c(0,4), xlim=c(-0.4,1.4), 
     main="College degree \n(0 = No, 1 = Yes)", lty=1, cex=1.2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==2,]$col, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==3,]$col, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==4,]$col, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE, family="Times")
legend("topleft", legend=c("Not Ally, Violates HR", "Not Ally, Protects HR", "Ally, Violates HR", "Ally, Protects HR"), lty=c(1,2,3,4), cex=1, bty="n")

# Age
plot(density(var_balance[var_balance$exp_4==1,]$midage, na.rm=T, bw=0.1), xlab="Age above 40 (0 = No, 1 = Yes)", ylim=c(0,4), xlim=c(-0.4,1.4), 
     main="Age \n(0 = Age<40, 1 = Age>=40)", lty=1, cex=1.2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==2,]$midage, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==3,]$midage, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==4,]$midage, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE, family="Times")
legend("topright", legend=c("Not Ally, Violates HR", "Not Ally, Protects HR", "Ally, Violates HR", "Ally, Protects HR"), lty=c(1,2,3,4), cex=1, bty="n")

# Party ID
plot(density(var_balance[var_balance$exp_4==1,]$republic, na.rm=T, bw=0.1), xlab="Republican (0 = No, 1 = Yes)", ylim=c(0,4), xlim=c(-0.4,1.4), 
     main="Party affiliation \n (0 = Non Republican, 1 = Republican)", lty=1, cex=1.2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==2,]$republic, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==3,]$republic, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==4,]$republic, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE, family="Times")
legend("topright", legend=c("Not Ally, Violates HR", "Not Ally, Protects HR", "Ally, Violates HR", "Ally, Protects HR"), lty=c(1,2,3,4), cex=1, bty="n")


# Democratic eval
plot(density(var_balance[var_balance$exp_4==1,]$dem_eval, na.rm=T, bw=0.1), xlab="Democratic evaluation >6 (0 = No, 1 = Yes)", ylim=c(0,4), xlim=c(-0.4,1.4), 
     main="Democratic evaluation \n(0 = Democ eval<6, 1=Democ eval>=6", lty=1, cex=1.2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==2,]$dem_eval, na.rm=T, bw=0.1), xlab="Democratic evaluation >6 (0 = No, 1 = Yes)", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==3,]$dem_eval, na.rm=T, bw=0.1), xlab="Democratic evaluation >6 (0 = No, 1 = Yes)", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(var_balance[var_balance$exp_4==4,]$dem_eval, na.rm=T, bw=0.1), xlab="Democratic evaluation >6 (0 = No, 1 = Yes)", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE, family="Times")
legend("topright", legend=c("Not Ally, Violates HR", "Not Ally, Protects HR", "Ally, Violates HR", "Ally, Protects HR"), lty=c(1,2,3,4), cex=1, bty="n")

# dev.off()

```

```{r, out.extra='angle=90', results='asis'}
#### Table S12: Regression Estimates of Support for War (Binary Dependent Variable, Pooled Sample) ####

# create a pooled sample
# load the yougov data
setwd("/Users/qingwang/Downloads/Data Replication")
df_ces <- read_dta("YouGov/YouGov_clean.dta")

# add dummy "purespec" and combine the two datasets
df_ps <- df_clean %>% 
  mutate(purespec = 1) %>%
  dplyr::select(attack, alliance, hmrts, male, age_cat, edu4, inc_10k, purespec, attack_cont)
df_ces <- df_ces %>% 
  mutate(purespec = 0) %>%
  dplyr::select(attack, alliance, hmrts, male, age_cat, edu4, inc_10k, purespec, attack_cont)
df_pool <- rbind(df_ps, df_ces)

# run the main analysis with binary DV
m1_pool <- lm(attack ~ hmrts + alliance + purespec, data = df_pool)
m2_pool <- lm(attack ~ alliance * hmrts + purespec, data = df_pool)
m3_pool <- lm(attack ~ alliance*hmrts + purespec +
           male + age_cat + edu4 + inc_10k, data = df_pool)

texreg(l = list(m1_pool, m2_pool, m3_pool),
       reorder.coef= c(2, 4, 5, 3, 6, 7, 8, 9, 1),
       custom.coef.names = c("(Intercept)", "Violating Human Rights", "PureSpectrum Sample",
                             "U.S. Military Alliance",
                             "Violating Human Rights $\\times$ U.S. Military Alliance",
                             "Male", "Age", "Education", "Income"),
       stars = c(0.01, 0.05, 0.1),
       digits = 4,
       caption = "Regression Estimates of Support for War (Binary Dependent Variable, Pooled Sample)",
       caption.above = T,
       include.ci = F,
       include.rmse =F,
       include.rsq = F,
       include.adjrs = F,
       label = "",
       custom.note = "",
       fontsize = "small") %>%
  gsub(".begin.center.", "\\\\centering", .) %>%
  gsub(".end.center.", "", .)

```

```{r}
#### Figure S11: Impact of Treatments on Support for War (Binary Dependent Variable, Pooled Sample, 95% Confidence Intervals) ####

### (S11a) Support for Attack

# create variable by treatment conditions
df_pool <- df_pool %>% 
  mutate(exp_4 = case_when(alliance == 0 & hmrts == 1 ~ 1,
                           alliance == 0 & hmrts == 0 ~ 2,
                           alliance == 1 & hmrts == 1 ~ 3,
                           alliance == 1 & hmrts == 0 ~ 4))

# calculate the mean support for attack in each treatment group
df_ate_pool <- df_pool %>% 
  group_by(exp_4, alliance, hmrts) %>%
  summarise(ate = mean(attack, na.rm = TRUE),
            n = n(),
            se = sd(attack, na.rm = TRUE) / sqrt(n)) %>%
  mutate(ci_low = ate - 1.96*se,
         ci_high = ate + 1.96*se)

# plot the mean support for attack in each treatment group and the 95% CI
p <- ggplot(df_ate_pool, aes(x = factor(exp_4, level=c('1', '3', '2', '4')), y = ate,
                        color=factor(exp_4))) + 
  theme_classic() +
  geom_point()+
  geom_errorbar(aes(ymin=ci_low, ymax=ci_high), width=.1,
                position=position_dodge(0.05)) +
  scale_x_discrete(labels= c('Non-Ally,Violates HR', 'Ally,Violates HR',
                             'Non-Ally,Protects HR', "Ally,Protects HR")) +
  scale_color_manual(values=c('#999999', 'black', 'black', 'black')) +
  theme(legend.position = "none") +
  labs(x = "", y = "Mean Support for Attack \n (% point)", size = 12) +
  geom_text(aes(label=round(ate, 1)), position=position_dodge(width=0.9), 
            vjust=.5, hjust = -.35) +
  theme(axis.text.x = element_text(angle = 20, hjust = 0.5, vjust = 0.5, size = 12),
        axis.title.y = element_text(size=12))

p
# ggsave("ate-pool.pdf", width = 6, height = 4)

### (S11b) Average Treatment Effect on Support for Attack
# calculate the difference in support for attack between 2-4 against baseline condition
est <- rep(NA, 4)
ci_low <- rep(NA, 4)
ci_high <- rep(NA, 4)
se <- rep(NA, 4)

for(i in 2:4){
  test <- t.test(df_pool$attack[df_pool$exp_4==i], 
                 df_pool$attack[df_pool$exp_4==1])
  est[i] <- test[["estimate"]][["mean of x"]] - test[["estimate"]][["mean of y"]]
  ci_low[i] <- test[["conf.int"]][1]
  ci_high[i] <- test[["conf.int"]][2]
  se[i] <- test[["stderr"]]
}

df_ate_diff_pool <- data.frame(exp_4 = df_ate_pool$exp_4, est, ci_low, ci_high, se)
df_ate_diff_pool <- df_ate_diff_pool[-1, ]

# plot the differences and 95% CI
p1 <- ggplot(df_ate_diff_pool, aes(x = factor(exp_4, level=c('3', '2', '4')), y = est)) + 
  theme_classic() +
  geom_point()+
  geom_errorbar(aes(ymin=ci_low, ymax=ci_high), width=.1,
                position=position_dodge(0.05)) +
  scale_x_discrete(labels= c('Ally,Violates HR', 'Non-Ally,Protects HR', 
                             'Ally,Protects HR')) +
  labs(x = "", y = "Average Treatment Effect \n (% point)", size = 12) +
  geom_text(aes(label=round(est, 1)), position=position_dodge(width=0.9), 
            vjust=.5, hjust = -.35) +
  geom_hline(yintercept = 0, linetype="dotted") +
  theme(axis.text.x = element_text(angle = 20, hjust = 0.5, vjust = 0.5, size = 12),
        axis.title.y = element_text(size=12))
p1
# ggsave("ate-diff-pool.pdf", width = 6, height = 4)

```

```{r, out.extra='angle=90', results='asis'}
#### Table S13: Regression Estimates of Support for War (5-Point Likert Scale, Pooled Sample) ####

# run the main analysis with binary DV
m1_cont_pool <- lm(attack_cont ~ hmrts + alliance + purespec, data = df_pool)
m2_cont_pool <- lm(attack_cont ~ alliance * hmrts + purespec, data = df_pool)
m3_cont_pool <- lm(attack_cont ~ alliance*hmrts + purespec +
                male + age_cat + edu4 + inc_10k, data = df_pool)

texreg(l = list(m1_cont_pool, m2_cont_pool, m3_cont_pool),
       reorder.coef= c(2, 4, 5, 3, 6, 7, 8, 9, 1),
       custom.coef.names = c("(Intercept)", "Violating Human Rights", "PureSpectrum Sample",
                             "U.S. Military Alliance",
                             "Violating Human Rights $\\times$ U.S. Military Alliance",
                             "Male", "Age", "Education", "Income"),
       stars = c(0.01, 0.05, 0.1),
       digits = 4,
       caption = "Regression Estimates of Support for War (5-point Likert Scale, Pooled Sample)",
       caption.above = T,
       include.ci = F,
       include.rmse =F,
       include.rsq = F,
       include.adjrs = F,
       label = "",
       custom.note = "",
       fontsize = "small") %>%
  gsub(".begin.center.", "\\\\centering", .) %>%
  gsub(".end.center.", "", .)

```


```{r}
#### Figure S12: Impact of Treatments on Support for War (5-Point Likert Scale, Pooled Sample, 95% Confidence Intervals) ####

### (S12a) Support for Attack

# calculate the mean support for attack in each treatment group
df_ate_cont_pool <- df_pool %>% 
  group_by(exp_4, alliance, hmrts) %>%
  summarise(ate = mean(attack_cont, na.rm = TRUE),
            n = n(),
            se = sd(attack, na.rm = TRUE) / sqrt(n)) %>%
  mutate(ci_low = ate - 1.96*se,
         ci_high = ate + 1.96*se)

# plot the mean support for attack in each treatment group and the 95% CI
p <- ggplot(df_ate_cont_pool, aes(x = factor(exp_4, level=c('1', '3', '2', '4')), y = ate,
                             color=factor(exp_4))) + 
  theme_classic() +
  geom_point()+
  geom_errorbar(aes(ymin=ci_low, ymax=ci_high), width=.1,
                position=position_dodge(0.05)) +
  scale_x_discrete(labels= c('Non-Ally,Violates HR', 'Ally,Violates HR',
                             'Non-Ally,Protects HR', "Ally,Protects HR")) +
  scale_color_manual(values=c('#999999', 'black', 'black', 'black')) +
  theme(legend.position = "none") +
  labs(x = "", y = "Mean Support for Attack \n (% point)", size = 12) +
  geom_text(aes(label=round(ate, 1)), position=position_dodge(width=0.9), 
            vjust=.5, hjust = -.35) +
  theme(axis.text.x = element_text(angle = 20, hjust = 0.5, vjust = 0.5, size = 12),
        axis.title.y = element_text(size=12))
p
# ggsave("ate_cont-pool.pdf", width = 6, height = 4)

### (S12b) Average Treatment Effect on Support for Attack
# calculate the difference in support for attack between 2-4 against baseline condition
est <- rep(NA, 4)
ci_low <- rep(NA, 4)
ci_high <- rep(NA, 4)
se <- rep(NA, 4)

for(i in 2:4){
  test <- t.test(df_pool$attack_cont[df_pool$exp_4==i], 
                 df_pool$attack_cont[df_pool$exp_4==1])
  est[i] <- test[["estimate"]][["mean of x"]] - test[["estimate"]][["mean of y"]]
  ci_low[i] <- test[["conf.int"]][1]
  ci_high[i] <- test[["conf.int"]][2]
  se[i] <- test[["stderr"]]
}

df_ate_cont_diff_pool <- data.frame(exp_4 = df_ate_cont_pool$exp_4, est, ci_low, ci_high, se)
df_ate_cont_diff_pool <- df_ate_cont_diff_pool[-1, ]


# plot the differences and 95% CI
p1 <- ggplot(df_ate_cont_diff_pool, aes(x = factor(exp_4, level=c('3', '2', '4')), y = est)) + 
  theme_classic() +
  geom_point()+
  geom_errorbar(aes(ymin=ci_low, ymax=ci_high), width=.1,
                position=position_dodge(0.05)) +
  scale_x_discrete(labels= c('Ally,Violates HR', 'Non-Ally,Protects HR', 
                             'Ally,Protects HR')) +
  labs(x = "", y = "Average Treatment Effect \n (% point)", size = 12) +
  geom_text(aes(label=round(est, 1)), position=position_dodge(width=0.9), 
            vjust=.5, hjust = -.35) +
  geom_hline(yintercept = 0, linetype="dotted") +
  theme(axis.text.x = element_text(angle = 20, hjust = 0.5, vjust = 0.5, size = 12),
        axis.title.y = element_text(size=12))

p1
# ggsave("ate_cont-diff-pool.pdf", width = 6, height = 4)

```


```{r,out.extra='angle=90', results='asis'}
#### Figure S13: The Estimated Effect of Alliance on Public Support for War (Pooled Sample) ####

# reproduce the estimates
f <- attack ~ hmrts*alliance
fit_pool  <- lm(f, data = df_pool)

# the smallest substantively interesting effect
pure_high_pool  <- 5.5
pure_low_pool  <- 11

# a data frame setting the values of the "other" variables
X_c <- data.frame(
  alliance = 0,  # low value
  hmrts = 1 # low value
)

quantile(df_pool$alliance, c(.1, .9))
# 10th percentile of alliance = 0
# 90th percentile of alliance = 1

# compute the comparison for alliance
alliance_comp_pool <- comparisons(fit_pool,
                             newdata = X_c, 
                             variables = list("alliance" = c(0, 1)), # low to high value
                             conf_level = 0.90)

# bind the comparisons together and plot
comp <- bind_rows(alliance_comp_pool , alliance_comp_pool)
comp <- comp %>% dplyr::mutate(ID = row_number(),
                               study = c("PureSpectrum Survey", 
                                         "PureSpectrum Survey"),
                               threshold = c("A Stringent Threshold", "A Baseline Threshold"),
                               m = c(pure_high_pool, pure_low_pool),
                               label = c("m = 5.5", "m = 11"))

set.seed(4755427)

# Create containers to store estimates and confidence intervals.
est.alliance_pool <- matrix(NA, nrow = 1, ncol = 3)

# Plot our estimate
est.alliance_pool[1, ] <- c(alliance_comp_pool$conf.low, alliance_comp_pool$estimate, alliance_comp_pool$conf.high)
rownames(est.alliance_pool) <- "The Estimate"

# pdf("equivalence-test-pool.pdf", 
#     height = 3, width = 9, family = "serif")
# par(mfrow = c(1, 1), oma = c(.5, 9, .5, 9), mar = c(1,1,1,1))
eplot(NULL, xlim = c(-12, 12), ylim = c(-1, 0), 
      xat = c(-11, -5.5, -1.9, 5.5, 11),
      anny = FALSE,
      xlab = "",
      xlabpos = 2.5)
abline(v = 11, xpd = FALSE, lty = 2, col = "blue")
abline(v = -11, xpd = FALSE, lty = 2, col = "blue")
abline(v = -5.5, xpd = FALSE, lty = 5, col = "red")
abline(v = 5.5, xpd = FALSE, lty = 5, col = "red")
legend(par('usr')[2], par('usr')[4], bty='n', xpd=NA,
       c("Stringent Threshold", "Baseline Threshold"), lty=c(5,2), col = c("red", "blue"))
for (i in 1:nrow(est.alliance_pool)) {
  est <- est.alliance_pool
  lines(c(est[i, 1], est[i, 3]), c(-i/(nrow(est) + 1), -i/(nrow(est) + 1)), lwd = 2) 
  points(est[i, 2], -i/(nrow(est) + 1), pch = 19, cex = .7)
  text(est[i, 2], -i/(nrow(est) + 1), rownames(est)[i], pos = 3, cex = .8)
}
# dev.off()

```

